\documentclass{article}
\usepackage[margin = 0.8in]{geometry}
\usepackage{amsmath,amssymb,amsthm,amsfonts}
\usepackage{graphicx}
% THEOREMS -------------------------------------------------------
\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition}
\theoremstyle{remark}
\newtheorem{rem}[thm]{Remark}
\numberwithin{equation}{section}
\newcommand{\var}{\operatorname{var}}
\newcommand{\mean}{\operatorname{mean}}
\newcommand{\cov}{\operatorname{cov}}
\newcommand{\dis}{\displaystyle}
% ----------------------------------------------------------------
\title{Control of surface roughness and RMS slope with patterned W and T}
\author{Jianqiao Huang}
\begin{document}
\maketitle
\tableofcontents
\section{Model Description}
In this work, we consider a 1D Edwards-Wilkinson (EW) type equation taking the following form
\begin{equation}
\frac{\partial h}{\partial t}=w+c_2\frac{\partial^2 h}{\partial x^2}+\xi(x,t) \label{eq:EW}
\end{equation}
where $x\in[0,L]$ is the spatial coordinates, $t$ is the time, $h(x,t)$ is the surface height, and $\xi(x,t)$ is a Gaussian white noise with a zero mean and the following covariance:
\begin{equation}
\left<\xi(x,t)\xi(x',t')\right>=\sigma^2\delta(x-x')\delta(t-t').
\end{equation}
where $\delta(\cdot)$ denotes the Dirac delta function. w is a patterned deposition profile with the form
\begin{equation}
\boxed{w(x) = w_0+A\sin\left(\frac{2k\pi}{L}x\right)}.
\end{equation}
The stochastic equation eq. \eqref{eq:EW} is subject to the following periodic boundary conditions
\begin{alignat}{2}
h(0,t) &= h(L,t)  \label{eq:EW_BC1}\\
\frac{\partial h}{\partial x}(0,t) &= \frac{\partial h}{\partial x}(L,t) \label{eq:EW_BC4}
\end{alignat}
and the initial condition
\begin{equation}
h(x,0)=h_0(x)
\end{equation}

To analyze the dynamics and obtain a finite-dimensional approximation of the EW equation, we first consider the eigenvalue problem of the linear operator of eq~\eqref{eq:EW} subject to the periodic boundary conditions of eq~\eqref{eq:EW_BC1}-\eqref{eq:EW_BC4}:
\begin{gather}
\mathcal{A}\phi_{n_x}(x) = c_2\frac{\partial^2\phi_{n_x}(x) }{\partial x^2}=\lambda_{n_x}\phi_{n_x}(x),\label{eq:EW_operator}\\
\nabla^j\phi_{n_x}(0)=\nabla^j\phi_{n_x}(L),\quad j = 0,1
\end{gather}
where $\lambda_{n_x}$ denotes an eigenvalue, $\phi_{n_x}$ denotes an eigenfunction, and $\nabla^j$, $j=0$, $1$, denotes the gradient of a given function. The solution of the eigenvalue problem is as follows:

\begin{align}
\lambda_{n_x}&= -\frac{4c_2\pi^2n_x^2}{L^2}\\
\phi_{1,n_x} =\phi_{n_x} &= \sqrt{\frac{2}{L}}\sin(\frac{2n_x\pi}{L}x)\\
\phi_{2,n_x} =\psi_{n_x} &= \left\{\begin{array}{ll}
                \sqrt{\frac{1}{L}}                               & \textrm{$n_x=0$ }\\
                \sqrt{\frac{2}{L}}\cos(\frac{2n_x\pi}{L}x)       & \textrm{$n_x\neq0$ }
                \end{array}\right.
\end{align}
The solution of the EW equation of eq~\eqref{eq:EW} can be expanded in an infinite series in terms of the eigenfunctions of the operator of eq~\eqref{eq:EW_operator} as follows:
\begin{equation}
h(x,t) = \sum_{n_x=0}^{+\infty}\phi_{1,n_x}z_{1,n_x}+\phi_{2,n_x}z_{2,n_x}=\sum_{n_x=0}^{+\infty}\phi_{n_x}\alpha_{n_x}+\psi_{n_x}\beta_{n_x},\label{eq:EW_solution}
\end{equation}
where $z_{1,n_x}(\alpha_{n_x})$, $z_{2,n_x}(\beta_{n_x})$ are time-varying coefficients.

Substituting the above expansion for the solution, $h(x,t)$, into eq~\eqref{eq:EW} and taking the inner product with the adjoint eigenfunctions, the following system of infinite stochastic linear ordinary differential equations (ODEs) for the temporal evolution of the time-varying coefficients is obtained:
\begin{gather}
\frac{dz_{2,0}}{dt} = \frac{d\beta_{0}}{dt}=w_{2,0}+\xi_{2,0}(t),\label{eq:mode0_ODE}\\
\frac{dz_{p,n_x}}{dt} = w_{p,n_x}+\lambda_{n_x}z_{p,n_x}+\xi_{p,n_x}(t)\\
p = 1,2, \quad n_x= 1,\cdots,\infty, \label{eq:mode_ODE}\notag
\end{gather}
where $\xi_{p,n_x}(t)=\int_0^{L}\xi(x,t)\phi_{p,n_x}(x)dx$ is the projection of the noise $\xi(x,t)$ on the ODE for $z_{p,n_x}$. The noise term, $\xi_{p,n_x}$, has zero mean and covariance
\begin{equation}
\left<\xi_{p,n_x}(t)\xi_{p,n_x}(t')\right>=\sigma^2\delta(t-t').
\end{equation}
Similarly, $w_{p,n_x}$ is the project of $w$ on the ODE for $z_{p,n_x}$, $w_{p,n_x}=\int_{0}^{L}\phi_{p,n_x}(x)w(x)dx$
\begin{itemize}
% w_{1,n_x} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{If $p=1$,
    \begin{itemize}
    \item{If $n_x\neq k$
            \begin{align}
            w_{1,n_x}     &=\boxed{0}
            \end{align}
        }

    \item{If $n_x= k$
        \begin{align}
        w_{1,n_x} &= A\sqrt{\frac{L}{2}}
        \end{align}
    }
    \end{itemize}
}
% w_{2,n_x} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{If $p=2$
\begin{itemize}
    \item{If $n_x\neq 0$
    \begin{align}
        w_{2,0}     &= \boxed{0}
    \end{align}
    }
    \item{If $n_x= 0$
    \begin{align}
    w_{2,0} &= w_0\sqrt{L}
    \end{align}
    }
\end{itemize}
}
\end{itemize}

The temporal evolution of the variance of mode $z_{p,n_x}$ can be obtained from the solution of the linear ODEs of eqs~\eqref{eq:mode0_ODE} and \eqref{eq:mode_ODE} as follows:
\begin{align}
\left<z_{2,0}(t)\right> &= w_{2,0}(t-t_0)\\
\var(z_{2,0}(t)        &= \sigma^2(t-t_0)\\
\left<z(t)\right>        &= e^{\lambda(t-t_0)}\left<z(t_0)\right>+\frac{w_p}{\lambda}(e^{\lambda(t-t_0)}-1)\\
\var(z(t))               &= e^{2\lambda(t-t_0)}\var(z(t_0))+\sigma^2\frac{e^{2\lambda(t-t_0)}-1}{2\lambda}\label{eq:cov_dyn}
\end{align}
where $z(t)=z_{p,n_x}(t)$ and $w_p=w_{p,n_x}$ for $n_x\neq0$.
% ====================================================================
\subsection{Surface roughness}
Surface roughness of the thin film is defined as the standard deviation of the surface height profile from its average height
\begin{equation}
r(t) = \sqrt{\frac{1}{L}\int_{0}^{L}\left[h(x,t)-\bar{h}(x,t)\right]^2dx}\label{eq:r_def}
\end{equation}
where $\bar{h}(t) = \frac{1}{L}\int_{0}^{L}h(x,t)dx$ is the average surface height. According to Eq.\eqref{eq:EW_solution}, we have
\begin{equation}
\bar{h}(t) = \frac{1}{L}\int_{0}^{L}\phi_{2,0}z_{2,0}dx=\sqrt{\frac{1}{L}}z_{2,0}.
\end{equation}
$\left<r^2(t)\right>$ can be rewritten as
\begin{align}
\left<r^2(t)\right>&=\left<\frac{1}{L}\int_{0}^{L}\left[h(x,t)-\bar{h}(x,t)\right]^2dx\right>\\
\end{align}
where
\begin{equation}
h(x,t)-\bar{h}(x,t)=\sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}{\phi_{p,n_x}z_{p,n_x}}
\end{equation}
\begin{align}
\begin{split}
\left<r^2(t)\right>&=\left<\frac{1}{L}\int_{0}^{L}\left[\sum_{p=1}^{2}\sum_{\substack{n_x=1}}^{\infty}z_{p,n_x}\phi_{p,n_x}(x)\right]^2dx\right>\\
                   &=\left<\frac{1}{L}\int_{0}^{L}\sum_{\substack{n_x=1}}^{\infty}\left(\phi_{1,n_x}^2z_{1,n_x}^2+\phi_{2,n_x}^2z_{2,n_x}^2\right)dx\right>\\
                   &=\frac{1}{L}\sum_{\substack{n_x=1}}^{\infty}\left(\left<z_{1,n_x}^2\right>+\left<z_{2,n_x}^2\right>\right)
\end{split}\label{eq:r2}
\end{align}
where
\begin{equation}
\boxed{\left<z^2_{p,n_x}\right>=\var(z_{p,n_x})+\left<z_{p,n_x}\right>^2\label{eq:mean_zp2}}.
\end{equation}
% ====================================================================
\subsection{RMS slope}
The rms slope is defined as the root--mean--square of the slope of the surface height
\begin{gather}
\begin{split}
m(t) &= \sqrt{\frac{1}{L}\int_0^{L}\left(\frac{\partial h}{\partial x}\right)^2dx}\\
     &= \sqrt{\frac{1}{L}\sum_{i=0}^{l-1}\left(\frac{\hat{h}(i+1)-\hat{h}(i)}{\Delta x}\right)^2\frac{L}{l}}
\end{split}\\
\begin{split}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\left<m^2(t)\right>&=\left<\frac{1}{L}\sum_{i=0}^{l-1}\left(\frac{\hat{h}(i+1)-\hat{h}(i)}{\Delta x}\right)^2\frac{L}{l}\right>\\
                   &=\left<\frac{1}{l\Delta x^2}\sum_{i=0}^{l-1}\left\{\sum_{p=1}^{2}\sum_{n_x=0}z_{p,n_x}\left[\phi_{p,n_x}(i+1)-\phi_{p,n_x}(i)\right]\right\}^2\right>\\
                   &=\left<\frac{1}{l\Delta x^2}\sum_{i=0}^{l-1}\sum_{p_1=1}^{2}\sum_{\substack{n_{1x}=0}}\sum_{p_2=1}^{2}\sum_{\substack{n_{2x}=0}}z_{p_1,n_{1x}}z_{p_2,n_{2x}}\Delta \phi_{p_1,n_{1x}}(i)\Delta\phi_{p_2,n_{2x}}(i)\right>\\
                   &=\frac{1}{l\Delta x^2}\sum_{p_1=1}^{2}\sum_{\substack{n_{1x}=0}}\sum_{p_2=1}^{2}\sum_{\substack{n_{2x}=0}}\left<z_{p_1,n_{1x}}z_{p_2,n_{2x}}\right>\left(\sum_{i=0}^{l-1}\Delta \phi_{p_1,n_{1x}}(i)\Delta\phi_{p_2,n_{2x}}(i)\right)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\left<m^2(t)\right>&=\left<\frac{1}{L^2}\iint_{0}^{L}\left[\sum_{p=1}^{4}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq0}}^{\infty}z_{p,n_x,n_y}\frac{\partial \phi_{p,n_x,n_y}(x,y)}{\partial x}\right]^2dxdy\right>\\
%                   &=\frac{1}{L^2}\iint_{0}^{L}\sum_{p_1=1}^{4}\sum_{p_2=1}^{4}\sum_{\substack{n_{1x},n_{1y}=0\\n_{1x}^2+n_{1y}^2\neq 0}}^{\infty}\sum_{\substack{n_{2x},n_{2y}=0\\n_{2x}^2+n_{2y}^2\neq 0}}^{\infty}\left<z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\right>\frac{\partial \phi_{p_1,n_{1x},n_{1y}}}{\partial x}\frac{\partial \phi_{p_2,n_{2x},n_{2y}}}{\partial x}dxdy\\
%                   &=\frac{1}{L^2}\sum_{p_1=1}^{4}\sum_{p_2=1}^{4}\sum_{\substack{n_{1x},n_{1y}=0\\n_{1x}^2+n_{1y}^2\neq 0}}^{\infty}\sum_{\substack{n_{2x},n_{2y}=0\\n_{2x}^2+n_{2y}^2\neq 0}}^{\infty}\left<z_{p_1,n_{1x},n_{1y}}z_{p_2,n_{2x},n_{2y}}\right>\iint_{0}^{L}\frac{\partial \phi_{p_1,n_{1x},n_{1y}}}{\partial x}\frac{\partial \phi_{p_2,n_{2x},n_{2y}}}{\partial x}dxdy\\
%                   &=\frac{1}{L^2}\sum_{p=1}^{4}\sum_{\substack{n_{x},n_{y}=0\\n_{x}^2+n_{y}^2\neq 0}}^{\infty}\left<z^2_{p,n_{x},n_{y}}\right>\iint_{0}^{L}\left(\frac{\partial \phi_{p,n_{x},n_{y}}}{\partial x}\right)^2dxdy\\
%                   &=\frac{1}{L^2}\sum_{\substack{n_x,n_y=0\\n_x^2+n_y^2\neq0}}^{\infty}\sum_{p=1}^{4}\left<z^2_{p,n_x,n_y}\right>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{split}\\
\begin{split}
 &\sum_{i=0}^{l-1}\Delta \phi_{p_1,n_{1x}}(i)\Delta\phi_{p_2,n_{2x}}(i)\\
=&\sum_{i=0}^{l-1}\left(\phi_{p_1,n_{1x}}(i+1)-\phi_{p_1,n_{1x}}(i)\right)\left(\phi_{p_2,n_{2x}}(i+1)-\phi_{p_2,n_{2x}}(i)\right)\\
=&\frac{2}{L}\left(\sum_{i=0}^{l-1}\left(\sin(\frac{2n_{1x}\pi}{l}(i+1))-\sin(\frac{2n_{1x}\pi}{l}i)\right)\left(\sin(\frac{2n_{2x}\pi}{l}(i+1))-\sin(\frac{2n_{2x}\pi}{l}i)\right)\right)\\
=&\frac{8}{L}\sin(\frac{n_1\pi}{l})\sin(\frac{n_2\pi}{l})\sum_{i=0}^{l-1}\left(\cos(\frac{n_1\pi}{l}(2i+1))\cos(\frac{n_2\pi}{l}(2i+1))\right)\\
\end{split}\\
\begin{split}
\left<m^2(t)\right> &=\frac{1}{l\Delta x^2}\sum_{p_1=1}^{2}\sum_{\substack{n_{1x}=0}}\sum_{p_2=1}^{2}\sum_{\substack{n_{2x}=0}}\left<z_{p_1,n_{1x}}z_{p_2,n_{2x}}\right>\left(\sum_{i=0}^{l-1}\Delta
                      \phi_{p_1,n_{1x}}(i)\Delta\phi_{p_2,n_{2x}}(i)\right)\\
                    &=\frac{1}{l\Delta x^2}\sum_{p=1}^{2}\sum_{\substack{n_x=0}}\left<z_{p,n_x}\right>^2\left(\frac{8}{L}\sin^2(\frac{n_x\pi}{l})\sum_{i=0}^{l-1}\left(\cos^2(\frac{n_x\pi}{l}(2i+1))\right)\right)
\end{split}\\
\left<m^2(t)\right> = \sum_{p=1}^{2}\sum_{\substack{n_{x}=0}}K_{p,n_x}\left<z_{p,n_x}^2\right>=\boxed{\sum_{p=1}^{2}\sum_{\substack{n_{x}=0}}\frac{4}{L\Delta x^2}\sin^2\left(\frac{\pi n_x}{l}\right)\left<z_{p,n_x}^2\right>}
\end{gather}
In terms of eigenfunctions
\begin{equation}
\left<m^2(t)\right>=\sum_{\substack{m=1}}^{l/2}\left(K_{1,m}\left<z^2_{1,m}\right>+K_{2,m}\left<z^2_{2,m}\right>\right)
\end{equation}
where
\begin{equation}
K_{p,m,n} = \
\end{equation}
% ==============================================================================
%\section{Detailed implementation}
%\section{Result validation}
\section{Fitting}
\subsection{Fitting to $r^2$ profile}
Equation~\eqref{eq:r2} and~\eqref{eq:mean_zp2} are used in the calculation of $\left<r^2\right>$ 
\begin{equation}
\begin{split}
\left<z^2_{p,n_x}\right> &= \var(z_{p,n_x})+\left<z_{p,n_x}\right>^2\\
                         &= 
\end{split}
\end{equation}

When $\lambda$ is small, 
\begin{equation}
\lim_{\lambda\rightarrow 0} \frac{e^{\lambda (t-t_0)}-1}{\lambda}=t-t_0 \label{eq:z_limit_0_lambda}
\end{equation}.
Thus 
\begin{equation}
\begin{split}
\left<z(t)\right>=e^{\lambda(t-t_0)}\left<z(t_0)\right>+w_p\left(t-t_0\right)
\end{split}
\end{equation}
\begin{equation}
\var(z(t))=e{2\lambda(t-t_0)}\var(z(t_0))+\sigma^2(t-t_0)
\end{equation}

Substitute eq.~\eqref{eq:z_limit_0_lambda} into eq~\eqref{eq:mean_zp2}
\begin{equation}

\end{equation}


% ==============================================================================
\section{MPC formulation and implementation using IPOPT}
MPC formulation similar to Ni Dong's and the one used in last IECR paper is used
\begin{gather}\label{eq:MPC}
\dis{\min_{T,W,A}} f(T,W,A) = q_{r^2}\left(r^2_{set}-\left<r^2_f\right>\right)^2+q_{m^2}\left(m^2_{set}-\left<m^2_f\right>\right)^2
\end{gather}
%\operatornamewithlimits{exp}_{x \in \mathbb{R}^n} J(t) = q_{r^2}\left(r^2(t_f)-r^2_{set}\right)^2+q_{m^2}\left(m^2(t_f)-m^2_{set}\right)^2\\
where
\begin{align}
\left<r^2_f\right> = \frac{1}{L}\sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}\left<z^2_{p,n_x}(t_f)\right>,\quad
\left<m^2_f\right> = \sum_{\substack{n_x=1}}\sum_{p=1}^{2}\left(K_{p,n_x}\left<z^2_{p,n_x}(t_f)\right>\right)
\end{align}
\begin{align}
\left<z^2_{p,n_x}(t_f)\right> &=\var(z_{p,n_x}(t_f))+\left<z_{p,n_x}(t_f)\right>^2\\
\left<z_{p,n_x}(t_f)\right>   &= e^{\lambda_{n_x}(t_f-t)}\left<z_{p,n_x}(t)\right>+\frac{w_p}{\lambda_{n_x}}(e^{\lambda_{n_x}(t_f-t)}-1)\\
\var(z_{p,n_x}(t_f))     &= e^{2\lambda_{n_x}(t_f-t)}\var(z_{p,n_x}(t))+\sigma^2\frac{e^{2\lambda_{n_x}(t_f-t)}-1}{2\lambda_{n_x}}\\
\lambda_{n_x}            &= -\frac{4c_2\pi^2}{L^2}n_x^2
\end{align}
\begin{align}
%c(T,W)        &= W\left(1-\frac{k_w}{W^{a_w}}e^{k_BT/E_w}\right)\\
c_2(W)      &= e^{aW^4+bW^3+cW^2+dW+e}\\
\sigma^2(W) &= fW+g
\end{align}
Subject to:
\begin{gather}
%T_{min}\leq T \leq T_{max}, \quad |T(t)-T(t-dt)|\leq \Delta T_{max}, \\
W_{min}\leq W \leq W_{max}, \quad |W(t)-W(t-dt)|\leq \Delta W_{max},\label{eq:MPC_LastConstraint}
\end{gather}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to implement in IPOPT, the optimization problem in eq.~\eqref{eq:MPC}--\eqref{eq:MPC_LastConstraint} is rewritten as
\begin{itemize}
\item{
Cost function
\begin{gather}\label{eq:MPC}
f(T,W,A) = q_{r^2}\left(r^2_{set}-\left<r^2_f\right>\right)^2+q_{m^2}\left(m^2_{set}-\left<m^2_f\right>\right)^2
\end{gather}
%\operatornamewithlimits{exp}_{x \in \mathbb{R}^n} J(t) = q_{r^2}\left(r^2(t_f)-r^2_{set}\right)^2+q_{m^2}\left(m^2(t_f)-m^2_{set}\right)^2\\
where
\begin{align}
\left<r^2_f\right> = \frac{1}{L}\sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}\left<z^2_{p,n_x}(t_f)\right>,\quad
\left<m^2_f\right> = \sum_{\substack{n_x=1}}\sum_{p=1}^{2}\left(K_{p,n_x}\left<z^2_{p,n_x}(t_f)\right>\right)
\end{align}
\begin{align}
\left<z^2_{p,n_x}(t_f)\right> &=\var(z_{p,n_x}(t_f))+\left<z_{p,n_x}(t_f)\right>^2\\
\left<z_{p,n_x}(t_f)\right>   &= e^{\lambda_{n_x}(t_f-t)}\left<z_{p,n_x}(t)\right>+\frac{w_p}{\lambda_{n_x}}(e^{\lambda_{n_x}(t_f-t)}-1)\\
\var(z_{p,n_x}(t_f))     &= e^{2\lambda_{n_x}(t_f-t)}\var(z_{p,n_x}(t))+\sigma^2\frac{e^{2\lambda_{n_x}(t_f-t)}-1}{2\lambda_{n_x}}\\
\lambda_{n_x}            &= -\frac{4c_2\pi^2}{L^2}n_x^2
\end{align}
\begin{align}
%c(T,W)        &= W\left(1-\frac{k_w}{W^{a_w}}e^{k_BT/E_w}\right)\\
c_2(W)      &= e^{aW^4+bW^3+cW^2+dW+e}\\
\sigma^2(W) &= fW+g
\end{align}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Gradient of cost function
\begin{gather}
\nabla f = \left[\frac{\partial f}{\partial W},\frac{\partial f}{\partial A}\right]
\end{gather}
%%%%%%%%%%%%%%%%%%% \frac{\partial f}{\partial T} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%$\frac{\partial f}{\partial T}$ is given by
%\begin{gather}
%\frac{\partial f}{\partial T} = -2q_{r^2}\left(r^2_{set}-r^2_f\right)\frac{\partial r^2_f}{\partial T}-2q_{m^2}\left(m^2_{set}-m^2_f\right)\frac{\partial m^2_f}{\partial T}
%\end{gather}
%\begin{gather}
%\frac{\partial r^2_f}{\partial T} = \frac{1}{L}\sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial T}\\
%\frac{\partial m^2_f}{\partial T} = \sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}\left(K_{p,n_x}\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial T}\right)
%\end{gather}
%\begin{gather}
%\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial T}=\frac{\partial \var(z_{p,n_x}(t_f))}{\partial T}+2\left<z_{p,n_x}\right>\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial T}
%\end{gather}
%\begin{gather}
%\frac{\partial \var(z_{p,n_x}(t_f))}{\partial T} = \frac{\partial \var(z_{p,n_x}(t_f))}{\partial \lambda_{n_x}}\frac{\partial \lambda_{n_x}}{\partial T}+\frac{\partial \var(z_{p,n_x}(t_f))}{\partial \sigma^2}\frac{\partial \sigma^2}{\partial T}
%\end{gather}
%\begin{gather}
%\begin{split}
%\frac{\partial \var(z_{p,n_x}(t_f))}{\partial \lambda_{n_x}} =& 2e^{2\lambda_{n_x}(t_f-t)}\var(z_{p,n_x}(t))(t_f-t)\\
%                                                                     &+\sigma^2\frac{4\lambda_{n_x}e^{2\lambda_{n_x}(t_f-t)}(t_f-t)-2\left(e^{2\lambda_{n_x}(t_f-t)}-1\right)}{4\lambda_{n_x}^2}
%\end{split}\\
%\frac{\partial \lambda_{n_x}}{\partial T}=-\frac{4\pi^2}{L^2}n_x^2\frac{\partial c_2}{\partial T}\\
%\frac{\partial c_2}{\partial T} = \frac{b}{l^2}e^{a+bT+cW}
%\end{gather}
%\begin{gather}
%\frac{\partial \var(z_{p,n_x}(t_f))}{\partial \sigma^2} = \frac{e^{2\lambda_{n_x}(t_f-t)}-1}{2\lambda_{n_x}}\\
%\frac{\partial \sigma^2}{\partial T} = d
%\end{gather}
%\begin{gather}
%\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial T}=\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial \lambda_{n_x}}\frac{\partial \lambda_{n_x}}{\partial T}\\
%\begin{split}
%\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial \lambda_{n_x}} &=e^{\lambda_{n_x}(t_f-t)}\left<z_{p,n_x}(t)\right>(t_f-t)\frac{\partial \lambda_{n_x}}{\partial T}\\
%&+w_p\frac{\lambda_{n_x}e^{\lambda_{n_x}(t_f-t)}(t_f-t)-(e^{\lambda_{n_x}(t_f-t)}-1)}{\lambda^2_{n_x}}
%\end{split}
%\end{gather}
%%%%%%%%%%%%%%%%%%% \frac{\partial f}{\partial W} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$\frac{\partial f}{\partial W}$ is given by
\begin{gather}
\frac{\partial f}{\partial W}= -2q_{r^2}\left(r^2_{set}-r^2_f\right)\frac{\partial r^2_f}{\partial W}-2q_{m^2}\left(m^2_{set}-m^2_f\right)\frac{\partial m^2_f}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial r^2_f}{\partial W} = \frac{1}{L}\sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial W}\\
\frac{\partial m^2_f}{\partial W} = \sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}\left(K_{p,n_x}\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial W}\right)
\end{gather}
\begin{gather}
\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial W}=\frac{\partial \var(z_{p,n_x}(t_f))}{\partial W}+2\left<z_{p,n_x}\right>\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial \var(z_{p,n_x}(t_f))}{\partial W} = \frac{\partial \var(z_{p,n_x}(t_f))}{\partial \lambda_{n_x}}\frac{\partial \lambda_{n_x}}{\partial W}+\frac{\partial \var(z_{p,n_x}(t_f))}{\partial \sigma^2}\frac{\partial \sigma^2}{\partial W}
\end{gather}
\begin{gather}
\frac{\partial \lambda_{n_x}}{\partial W}= -\frac{4\pi^2}{L^2}n_x^2\frac{\partial c_2}{\partial W}\\
\frac{\partial c_2}{\partial W} = \left(4aW^3+3bW^2+2cW+d\right)e^{aW^4+bW^3+cW^2+dW+e}
\end{gather}
\begin{gather}
\frac{\partial \sigma^2}{\partial W} = f
\end{gather}
\begin{gather}
\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial W}=\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial \lambda_{n_x}}\frac{\partial \lambda_{n_x}}{\partial W}+\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial w_{p,m}}\frac{\partial w_{p,m}}{\partial W}\\
\frac{\partial w_{p,m}}{\partial W}=
\begin{cases}
L    & p=2,m=0\\
0    & \text{otherwise}
\end{cases}
\end{gather}
%%%%%%%%%%%%%%%%%%% \frac{\partial f}{\partial A} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$\frac{\partial f}{\partial A}$ is given by
\begin{gather}
\frac{\partial f}{\partial A}= -2q_{r^2}\left(r^2_{set}-r^2_f\right)\frac{\partial r^2_f}{\partial A}-2q_{m^2}\left(m^2_{set}-m^2_f\right)\frac{\partial m^2_f}{\partial A}
\end{gather}
\begin{gather}
\frac{\partial r^2_f}{\partial A} = \frac{1}{L}\sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial A}\\
\frac{\partial m^2_f}{\partial A} = \sum_{\substack{n_x=1}}^{\infty}\sum_{p=1}^{2}\left(K_{p,n_x}\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial A}\right)
\end{gather}
\begin{gather}
\frac{\partial \left<z^2_{p,n_x}(t_f)\right>}{\partial A}=2\left<z_{p,n_x}\right>\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial A}
\end{gather}
\begin{gather}
\begin{split}
\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial A}=\frac{\partial \left<z_{p,n_x}(t_f)\right>}{\partial w_p}\frac{\partial w_p}{\partial A} = \frac{1}{\lambda_{n_x}}\left(e^{\lambda_{n_x}(t_f-t)}-1\right)\frac{\partial w_{p,n_x}}{\partial A}
\end{split}\\
\frac{\partial w_{p,n_x}}{\partial A} =
\begin{cases}
0                                         & \text{if $n_x\neq k$}\\
\sqrt{\frac{L}{2}} & \text{if $p=1, n_x=k$}\\
\end{cases}
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Inequality constraints
\begin{gather}
0\leq g(x) = W-A \leq \infty
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Jacobian matrix of constraints
\begin{gather}
\nabla g = \begin{bmatrix} 1 & -1 \end{bmatrix}
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Variable constraints
\begin{gather}
%\max\{T_{min},T_{old}-\Delta T_{max}\}\leq T \leq \min\{T_{max},T_{old}+\Delta T_{max}\}\\
\max\{W_{min},W_{old}-\Delta W_{max}\}\leq W \leq \min\{W_{max},W_{old}+\Delta W_{max}\}\\
0 \leq A \leq \infty
\end{gather}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\item{
Hessian matrix of the Lagrangian
\begin{gather}
\sigma_f\nabla^2 f+\lambda_1\nabla^2 g=
\sigma_f
\begin{bmatrix}
%\frac{\partial^2 f}{\partial T^2}         &                                            & \\
%\frac{\partial^2 f}{\partial W\partial T} & \frac{\partial^2 f}{\partial W^2}          & \\
%\frac{\partial^2 f}{\partial A\partial T} & \frac{\partial^2 f}{\partial A\partial T}  & \frac{\partial^2 f}{\partial A^2}
\frac{\partial^2 f}{\partial W^2 } & \\
\frac{\partial^2 f}{\partial W\partial A}  & \frac{\partial^2 f}{\partial A^2}
\end{bmatrix}
\end{gather}
The Hessian matrix will be approximated by finite difference method (quasi Newton method), thus there is no need for detail analytical expression. The Hessian matrix has 3 non-zero elements.
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{itemize}
\end{document}
